home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zgebrd.f < prev    next >
Text File  |  1996-07-19  |  9KB  |  260 lines
  1.       SUBROUTINE ZGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
  2.      $                   INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            INFO, LDA, LWORK, M, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   D( * ), E( * )
  14.       COMPLEX*16         A( LDA, * ), TAUP( * ), TAUQ( * ),
  15.      $                   WORK( LWORK )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  ZGEBRD reduces a general complex M-by-N matrix A to upper or lower
  22. *  bidiagonal form B by a unitary transformation: Q**H * A * P = B.
  23. *
  24. *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
  25. *
  26. *  Arguments
  27. *  =========
  28. *
  29. *  M       (input) INTEGER
  30. *          The number of rows in the matrix A.  M >= 0.
  31. *
  32. *  N       (input) INTEGER
  33. *          The number of columns in the matrix A.  N >= 0.
  34. *
  35. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  36. *          On entry, the M-by-N general matrix to be reduced.
  37. *          On exit,
  38. *          if m >= n, the diagonal and the first superdiagonal are
  39. *            overwritten with the upper bidiagonal matrix B; the
  40. *            elements below the diagonal, with the array TAUQ, represent
  41. *            the unitary matrix Q as a product of elementary
  42. *            reflectors, and the elements above the first superdiagonal,
  43. *            with the array TAUP, represent the unitary matrix P as
  44. *            a product of elementary reflectors;
  45. *          if m < n, the diagonal and the first subdiagonal are
  46. *            overwritten with the lower bidiagonal matrix B; the
  47. *            elements below the first subdiagonal, with the array TAUQ,
  48. *            represent the unitary matrix Q as a product of
  49. *            elementary reflectors, and the elements above the diagonal,
  50. *            with the array TAUP, represent the unitary matrix P as
  51. *            a product of elementary reflectors.
  52. *          See Further Details.
  53. *
  54. *  LDA     (input) INTEGER
  55. *          The leading dimension of the array A.  LDA >= max(1,M).
  56. *
  57. *  D       (output) DOUBLE PRECISION array, dimension (min(M,N))
  58. *          The diagonal elements of the bidiagonal matrix B:
  59. *          D(i) = A(i,i).
  60. *
  61. *  E       (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
  62. *          The off-diagonal elements of the bidiagonal matrix B:
  63. *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
  64. *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
  65. *
  66. *  TAUQ    (output) COMPLEX*16 array dimension (min(M,N))
  67. *          The scalar factors of the elementary reflectors which
  68. *          represent the unitary matrix Q. See Further Details.
  69. *
  70. *  TAUP    (output) COMPLEX*16 array, dimension (min(M,N))
  71. *          The scalar factors of the elementary reflectors which
  72. *          represent the unitary matrix P. See Further Details.
  73. *
  74. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  75. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  76. *
  77. *  LWORK   (input) INTEGER
  78. *          The length of the array WORK.  LWORK >= max(1,M,N).
  79. *          For optimum performance LWORK >= (M+N)*NB, where NB
  80. *          is the optimal blocksize.
  81. *
  82. *  INFO    (output) INTEGER
  83. *          = 0:  successful exit.
  84. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  85. *
  86. *  Further Details
  87. *  ===============
  88. *
  89. *  The matrices Q and P are represented as products of elementary
  90. *  reflectors:
  91. *
  92. *  If m >= n,
  93. *
  94. *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
  95. *
  96. *  Each H(i) and G(i) has the form:
  97. *
  98. *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
  99. *
  100. *  where tauq and taup are complex scalars, and v and u are complex
  101. *  vectors; v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in
  102. *  A(i+1:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in
  103. *  A(i,i+2:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
  104. *
  105. *  If m < n,
  106. *
  107. *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
  108. *
  109. *  Each H(i) and G(i) has the form:
  110. *
  111. *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
  112. *
  113. *  where tauq and taup are complex scalars, and v and u are complex
  114. *  vectors; v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in
  115. *  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in
  116. *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
  117. *
  118. *  The contents of A on exit are illustrated by the following examples:
  119. *
  120. *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
  121. *
  122. *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
  123. *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
  124. *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
  125. *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
  126. *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
  127. *    (  v1  v2  v3  v4  v5 )
  128. *
  129. *  where d and e denote diagonal and off-diagonal elements of B, vi
  130. *  denotes an element of the vector defining H(i), and ui an element of
  131. *  the vector defining G(i).
  132. *
  133. *  =====================================================================
  134. *
  135. *     .. Parameters ..
  136.       COMPLEX*16         ONE
  137.       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
  138. *     ..
  139. *     .. Local Scalars ..
  140.       INTEGER            I, IINFO, J, LDWRKX, LDWRKY, MINMN, NB, NBMIN,
  141.      $                   NX
  142.       DOUBLE PRECISION   WS
  143. *     ..
  144. *     .. External Subroutines ..
  145.       EXTERNAL           XERBLA, ZGEBD2, ZGEMM, ZLABRD
  146. *     ..
  147. *     .. Intrinsic Functions ..
  148.       INTRINSIC          MAX, MIN
  149. *     ..
  150. *     .. External Functions ..
  151.       INTEGER            ILAENV
  152.       EXTERNAL           ILAENV
  153. *     ..
  154. *     .. Executable Statements ..
  155. *
  156. *     Test the input parameters
  157. *
  158.       INFO = 0
  159.       IF( M.LT.0 ) THEN
  160.          INFO = -1
  161.       ELSE IF( N.LT.0 ) THEN
  162.          INFO = -2
  163.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  164.          INFO = -4
  165.       ELSE IF( LWORK.LT.MAX( 1, M, N ) ) THEN
  166.          INFO = -10
  167.       END IF
  168.       IF( INFO.LT.0 ) THEN
  169.          CALL XERBLA( 'ZGEBRD', -INFO )
  170.          RETURN
  171.       END IF
  172. *
  173. *     Quick return if possible
  174. *
  175.       MINMN = MIN( M, N )
  176.       IF( MINMN.EQ.0 ) THEN
  177.          WORK( 1 ) = 1
  178.          RETURN
  179.       END IF
  180. *
  181.       WS = MAX( M, N )
  182.       LDWRKX = M
  183.       LDWRKY = N
  184. *
  185. *     Set the block size NB and the crossover point NX.
  186. *
  187.       NB = MAX( 1, ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 ) )
  188. *
  189.       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
  190. *
  191. *        Determine when to switch from blocked to unblocked code.
  192. *
  193.          NX = MAX( NB, ILAENV( 3, 'ZGEBRD', ' ', M, N, -1, -1 ) )
  194.          IF( NX.LT.MINMN ) THEN
  195.             WS = ( M+N )*NB
  196.             IF( LWORK.LT.WS ) THEN
  197. *
  198. *              Not enough work space for the optimal NB, consider using
  199. *              a smaller block size.
  200. *
  201.                NBMIN = ILAENV( 2, 'ZGEBRD', ' ', M, N, -1, -1 )
  202.                IF( LWORK.GE.( M+N )*NBMIN ) THEN
  203.                   NB = LWORK / ( M+N )
  204.                ELSE
  205.                   NB = 1
  206.                   NX = MINMN
  207.                END IF
  208.             END IF
  209.          END IF
  210.       ELSE
  211.          NX = MINMN
  212.       END IF
  213. *
  214.       DO 30 I = 1, MINMN - NX, NB
  215. *
  216. *        Reduce rows and columns i:i+ib-1 to bidiagonal form and return
  217. *        the matrices X and Y which are needed to update the unreduced
  218. *        part of the matrix
  219. *
  220.          CALL ZLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
  221.      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
  222.      $                WORK( LDWRKX*NB+1 ), LDWRKY )
  223. *
  224. *        Update the trailing submatrix A(i+ib:m,i+ib:n), using
  225. *        an update of the form  A := A - V*Y' - X*U'
  226. *
  227.          CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-I-NB+1,
  228.      $               N-I-NB+1, NB, -ONE, A( I+NB, I ), LDA,
  229.      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
  230.      $               A( I+NB, I+NB ), LDA )
  231.          CALL ZGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
  232.      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
  233.      $               ONE, A( I+NB, I+NB ), LDA )
  234. *
  235. *        Copy diagonal and off-diagonal elements of B back into A
  236. *
  237.          IF( M.GE.N ) THEN
  238.             DO 10 J = I, I + NB - 1
  239.                A( J, J ) = D( J )
  240.                A( J, J+1 ) = E( J )
  241.    10       CONTINUE
  242.          ELSE
  243.             DO 20 J = I, I + NB - 1
  244.                A( J, J ) = D( J )
  245.                A( J+1, J ) = E( J )
  246.    20       CONTINUE
  247.          END IF
  248.    30 CONTINUE
  249. *
  250. *     Use unblocked code to reduce the remainder of the matrix
  251. *
  252.       CALL ZGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
  253.      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
  254.       WORK( 1 ) = WS
  255.       RETURN
  256. *
  257. *     End of ZGEBRD
  258. *
  259.       END
  260.